%% 计算每个区域内的数据点的q值
function [Q]=area_Q(DATA,HEX,track_hex,track_locat,lat1,lat2,lon1,lon2,xx) %步长a; %经纬度范围
% DATA = evalin('base','DATA');  %提取DATA
% HEX = evalin('base','HEX');  %提取HEX

q=[0]; %创建一个矩阵q存储每个小区域中的Q值：
m=1;n=0;t=1;k=[0];
for lat=lat1:xx:lat2-xx %步长0.5分区域
    n=0;
    for lon=lon1:xx:lon2-xx
        k(t)=1;
        for i=1:1:length(DATA) %遍历数据点 
             if DATA(i,2)>=lat && DATA(i,2)<lat+xx && DATA(i,1)>=lon && DATA(i,1)<lon+xx
                 [v_next]=next_position1(i,DATA,HEX); 
                %计算三个参数：v;nic;
                v=1-1/exp((abs(DATA(i,4)-v_next))/20);
                nic=1-DATA(i,3)/9;
                q(t,k(t))=0.2*v+0.8*nic; %计算对应的q值
                k(t)=k(t)+1;                  
             end
        end 
        n=n+1;  
        t=t+1;       
    end 
    m=m+1; 
end
%% 轨迹断裂处补点
% track_hex = evalin('base','track_hex');  %提取异常轨迹数据
% track_locat = evalin('base','track_locat');  

DATA_1=[];
[x,y]=line_point(track_locat(1,2:3),track_locat(1,4:5));
DATA_1(:,1)=x';DATA_1(:,2)=y';
 for i=2:1:length(track_locat)
     [x,y]=line_point(track_locat(i,2:3),track_locat(i,4:5));
     t=[x',y'];
     DATA_1=[DATA_1;t];
 end
%% 补点的q值
m=1;n=0;t=1;
for lat=lat1:xx:lat2-xx
    n=0;
    for lon=lon1:xx:lon2-xx
        for i=1:1:length(DATA_1)   
             if DATA_1(i,2)>=lat && DATA_1(i,2)<lat+xx && DATA_1(i,1)>=lon && DATA_1(i,1)<lon+xx
                q(t,k(t))=1; %对应的q=1
                k(t)=k(t)+1;                  
             end
        end 
        n=n+1; 
        t=t+1;       
    end 
    m=m+1; 
end
%% 计算所有点的Q值和中误差
Q=[0]; Q_1=[0];%存储Q值,中误差值
m=1;n=0;t=1;
for lat=lat1:xx:lat2-xx %步长0.5分区域
    n=0;
    for lon=lon1:xx:lon2-xx
        n=n+1;  
        if k(t)>=10 %区域内数据点数>10
            Q(m,n)=mean(q(t,1:k(t)-1));  %计算Q
            %计算中误差Q_1
            q_2=0;
            for a=1:1:k(t)-1
                qq=(q(t,a)-Q(m,n))^2;
                q_2=q_2+qq;
            end
            Q_1(m,n)=sqrt(q_2)/(k(t)-1);  %Q_1为中误差 
        else %区域内数据点数<10
            Q(m,n)=-1;
            Q_1(m,n)=-1;
        end 
        t=t+1;       
    end 
    m=m+1; 
end
%% 检测区域可视化
m=1;n=0;t=1;
figure(1);
for lat=lat1:xx:lat2-xx %步长0.5分区域
    n=0;
    for lon=lon1:xx:lon2-xx
        n=n+1; 
        if Q(m,n)>=0.9 %Q>0.9;红色填充区域
            X = [lon,lon,lon+xx,lon+xx];
            Y = [lat,lat+xx,lat+xx,lat];
            fill(X,Y,'r')
            hold on
        end
        if Q(m,n)>=0.5 && Q(m,n)<0.9 %Q>0.5 & Q<0.9;黄色填充区域
            X = [lon,lon,lon+xx,lon+xx];
            Y = [lat,lat+xx,lat+xx,lat];
            fill(X,Y,'y')
            hold on
        end
        if Q(m,n)<0.5 && Q(m,n)>=0 %Q<0.5 & Q>0.9;绿色填充区域
            X = [lon,lon,lon+xx,lon+xx];
            Y = [lat,lat+xx,lat+xx,lat];
            fill(X,Y,'g')
            hold on
        end
        if Q(m,n)<0 %Q<0 ;蓝色填充区域
            X = [lon,lon,lon+xx,lon+xx];
            Y = [lat,lat+xx,lat+xx,lat];
            fill(X,Y,'b')
            hold on
        end
        t=t+1;       
    end 
    m=m+1; 
end
hold on
%绘制对应区域的数据点分布图
% plot(DATA(:,1), DATA(:,2), '.b');
grid on
hold on
title('区域检测结果')
xlim([lon1 lon2]);      %右半边坐标轴的显示范围
ylim([lat1 lat2]);       %右半边坐标轴的显示范围
xlabel('经度')
ylabel('纬度')

end


